function llhess = HessPhat(theta,~,Y,X,Z,YY,V_R,V_US,V_UK,V_France,V_Russia,V_China,U0,intervEq,DMES,reb,mpow)

N=size(X,1);
KZ=size(Z,2)/5;
NS=size(V_R,2);
nBeta=length(reb{1});
nPhi=length(mpow{1});
nGamma=length(reb{2});
nChi=length(mpow{2});
nYY=size(YY,1);

Xrebel=X(:,reb{1}); Xmpow=X(:,mpow{1});

ZUS=Z(:,1:KZ);
ZUK=Z(:,KZ+(1:KZ)); 
ZFRN=Z(:,2*KZ+(1:KZ));
ZRUS=Z(:,3*KZ+(1:KZ));
ZCHN=Z(:,4*KZ+(1:KZ));

ZrebelUS=ZUS(:,reb{2}); ZmpowUS=ZUS(:,mpow{2});
ZrebelUK=ZUK(:,reb{2}); ZmpowUK=ZUK(:,mpow{2});
ZrebelFRN=ZFRN(:,reb{2}); ZmpowFRN=ZFRN(:,mpow{2});
ZrebelRUS=ZRUS(:,reb{2}); ZmpowRUS=ZRUS(:,mpow{2});
ZrebelCHN=ZCHN(:,reb{2}); ZmpowCHN=ZCHN(:,mpow{2});

beta=theta(1:nBeta,1);

phi_US=theta(nBeta+[1,5+(1:nPhi-1)],1);
phi_UK=theta(nBeta+[2,5+(1:nPhi-1)],1);
phi_France=theta(nBeta+[3,5+(1:nPhi-1)],1);
phi_Russia=theta(nBeta+[4,5+(1:nPhi-1)],1);
phi_China=theta(nBeta+[5,5+(1:nPhi-1)],1);
    
gamma_US=theta(nBeta+4+nPhi+1,1);
gamma_UK=theta(nBeta+4+nPhi+2,1);
gamma_France=theta(nBeta+4+nPhi+3,1);
gamma_Russia=theta(nBeta+4+nPhi+4,1);
gamma_China=theta(nBeta+4+nPhi+5,1);
gamma0=theta(nBeta+4+nPhi+5+(1:nGamma),1);

chi=theta(nBeta+4+nPhi+5+nGamma+(1:nChi),1);

delta_US_UK=theta(nBeta+4+nPhi+5+nGamma+nChi+1,1);
delta_US_France=theta(nBeta+4+nPhi+5+nGamma+nChi+2,1);
delta_US_Russia=theta(nBeta+4+nPhi+5+nGamma+nChi+3,1);
delta_US_China=theta(nBeta+4+nPhi+5+nGamma+nChi+4,1);
delta_UK_France=theta(nBeta+4+nPhi+5+nGamma+nChi+5,1);
delta_UK_Russia=theta(nBeta+4+nPhi+5+nGamma+nChi+6,1);
delta_UK_China=theta(nBeta+4+nPhi+5+nGamma+nChi+7,1);
delta_France_Russia=theta(nBeta+4+nPhi+5+nGamma+nChi+8,1);
delta_France_China=theta(nBeta+4+nPhi+5+nGamma+nChi+9,1);
delta_Russia_China=theta(nBeta+4+nPhi+5+nGamma+nChi+10,1);

lambda=theta(nBeta+4+nPhi+5+nGamma+nChi+11:end,1);

llhess=zeros(length(theta));

for n=1:N
    
    u_R = YY(:,1) .* (Xrebel(n,:)*beta+YY(:,2:end)*[gamma_US+ZrebelUS(n,:)*gamma0;gamma_UK+ZrebelUK(n,:)*gamma0;...
        gamma_France+ZrebelFRN(n,:)*gamma0;gamma_Russia+ZrebelRUS(n,:)*gamma0;gamma_China+ZrebelCHN(n,:)*gamma0]);
    
    u_US = YY(:,2) .* (Xmpow(n,:)*phi_US + ZmpowUS(n,:)*chi +...
        YY(:,3:6) * [delta_US_UK; delta_US_France; delta_US_Russia; delta_US_China]);
    
    u_UK = YY(:,3) .* (Xmpow(n,:)*phi_UK + ZmpowUK(n,:)*chi +...
        YY(:,[2,4:6])*[delta_US_UK; delta_UK_France; delta_UK_Russia; delta_UK_China]);
    
    u_France = YY(:,4) .* (Xmpow(n,:)*phi_France + ZmpowFRN(n,:)*chi +...
        YY(:,[2:3,5:6])*[delta_US_France; delta_UK_France; delta_France_Russia; delta_France_China]);
    
    u_Russia = YY(:,5) .* (Xmpow(n,:)*phi_Russia + ZmpowRUS(n,:)*chi +...
        YY(:,[2:4,6])*[delta_US_Russia; delta_UK_Russia; delta_France_Russia; delta_Russia_China]);
    
    u_China = YY(:,6) .* (Xmpow(n,:)*phi_China + ZmpowCHN(n,:)*chi +...
        YY(:,2:5)*[delta_US_China; delta_UK_China; delta_France_China; delta_Russia_China]);
       
    V=[V_R(:,:,n);V_US(:,:,n);V_UK(:,:,n);V_France(:,:,n);V_Russia(:,:,n);V_China(:,:,n)]';
    U=[u_R;u_US;u_UK;u_France;u_Russia;u_China]';
    U0n=reshape(U0(:,:,n),6*nYY,1)';
    weights = exp(-sum((V-U).^2,2)/2 + sum((V-U0n).^2,2)/2) / NS;
    
    %Gradients w.r.t. payoffs
    Du=zeros(6*nYY,length(theta)-length(lambda));
    % rebel payoffs  
    Du(1:nYY,[1:nBeta,nBeta+4+nPhi+(1:5+nGamma)])=[YY(:,1)*Xrebel(n,:),YY(:,2:end),YY(:,2)*ZrebelUS(n,:)+...
    	YY(:,3)*ZrebelUK(n,:)+YY(:,4)*ZrebelFRN(n,:)+YY(:,5)*ZrebelRUS(n,:)+YY(:,6)*ZrebelCHN(n,:)];
	% US payoffs
    Du(nYY+(1:nYY),[nBeta+[1,5+(1:nPhi-1)],nBeta+4+nPhi+5+nGamma+(1:nChi),nBeta+4+nPhi+5+nGamma+nChi+(1:4)])=YY(:,2).*[YY(:,1)*[Xmpow(n,:),ZmpowUS(n,:)],YY(:,3:6)];
    % UK
    Du(2*nYY+(1:nYY),[nBeta+[2,5+(1:nPhi-1)],nBeta+4+nPhi+5+nGamma+(1:nChi),nBeta+4+nPhi+5+nGamma+nChi+[1,5:7]])=YY(:,3).*[YY(:,1)*[Xmpow(n,:),ZmpowUK(n,:)],YY(:,[2,4:6])];
    % France 
    Du(3*nYY+(1:nYY),[nBeta+[3,5+(1:nPhi-1)],nBeta+4+nPhi+5+nGamma+(1:nChi),nBeta+4+nPhi+5+nGamma+nChi+[2,5,8:9]])=YY(:,4).*[YY(:,1)*[Xmpow(n,:),ZmpowFRN(n,:)],YY(:,[2:3,5:6])];
    % Russia
    Du(4*nYY+(1:nYY),[nBeta+[4,5+(1:nPhi-1)],nBeta+4+nPhi+5+nGamma+(1:nChi),nBeta+4+nPhi+5+nGamma+nChi+[3,6,8,10]])=YY(:,5).*[YY(:,1)*[Xmpow(n,:),ZmpowRUS(n,:)],YY(:,[2:4,6])];
    % China
    Du(5*nYY+(1:nYY),[nBeta+[5,5+(1:nPhi-1)],nBeta+4+nPhi+5+nGamma+(1:nChi),nBeta+4+nPhi+5+nGamma+nChi+[4,7,9:10]])=YY(:,6).*[YY(:,1)*[Xmpow(n,:),ZmpowCHN(n,:)],YY(:,2:5)];

    Pra = zeros(NS, 1);
    Grad = zeros(NS, length(lambda));
    H22=zeros(length(lambda));
	for b=1:NS
    	[Pra(b,1), Grad(b,:), h22] = profprob(Y(n,1), DMES{b,n}, intervEq{b,n}, V_R(1,b,n), lambda);
        H22 = H22 + weights(b,1) * h22;
	end
    
    ln = weights' * Pra;
    lg1 = weights' * (Pra .* (V-U) * Du) / ln;
    lg2 = weights' * Grad / ln;
    
    H11 = - Du' * Du + (weights .* (V-U) * Du)' * (Pra .* (V-U) * Du) / ln - lg1' * lg1;
    H12 = (weights .* (V-U) * Du)' * Grad / ln - lg1' * lg2;
    H22 = H22 / ln - lg2' * lg2;
    
    llhess = llhess - [H11, H12; H12', H22];
        
end         
    


